Liquid argon benchmarks

Sebastian Micluța-Câmpeanu, Mikhail Vaganov
using ProgressLogging
using NBodySimulator, OrdinaryDiffEq, StaticArrays
using Plots, DataFrames, StatsPlots

function setup(t)
    T = 120.0 # K
    kb = 1.38e-23 # J/K
    ϵ = T * kb # J
    σ = 3.4e-10 # m
    ρ = 1374 # kg/m^3
    m = 39.95 * 1.6747 * 1e-27 # kg
    N = 216
    L = (m*N/ρ)^(1/3)
    R = 2.25σ
    v_dev = sqrt(kb * T / m) # m/s

    _L = L / σ
     = 1.0
     = 1.0
    _m = 1.0
    _v = v_dev / sqrt(ϵ / m)
    _R = R / σ

    bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
    lj_parameters = LennardJonesParameters(, , _R)
    pbc = CubicPeriodicBoundaryConditions(_L)
    lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters));
    simulation = NBodySimulation(lj_system, (0.0, t), pbc, /T)

    return simulation
end
setup (generic function with 1 method)

We will consider the following integrators

config(τ, at, rt) = [
    # symplectic
    (alg=VelocityVerlet, dt=τ),
    (alg=McAte2, dt=τ),
    # (alg=CalvoSanz4, dt=τ),
    # (alg=McAte5, dt=τ),
    # (alg=Yoshida6, dt=τ),
    # (alg=KahanLi8, dt=τ),
    # DPRKN
    (alg=DPRKN6, abstol=at, rtol=rt),
    # (alg=DPRKN8, abstol=at, rtol=rt),
    (alg=DPRKN12, abstol=at, rtol=rt),
    # others
    (alg=Tsit5, abstol=at, rtol=rt),
    (alg=Vern7, abstol=at, rtol=rt),
    (alg=Vern9, abstol=at, rtol=rt)
]
config (generic function with 1 method)

We will first begin with the evolution of the energy error over time for different integrators. In this case we will fix both simulation time and timestep.

function benchmark(energyerr, rts, ts, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE(t) = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = [ΔE(t) for t in sol.t[2:end]]
        rts[alg] = rt
        ts[alg] = sol.t[2:end]
    end
end

ΔE = Dict()
rt = Dict()
ts = Dict()
configs = config(2.3e-4, 1e-20, 1e-20)
benchmark(ΔE, rt, ts, 40., configs)

plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:topleft);
for c in configs
    plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", xscale=:log10, yscale=:log10)
end
plt

Now, let us consider a fixed simulation time

t = 35.0
results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[], :abstol=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[]);

function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt, b, gc, memalloc = @timed solve(prob, alg();
            save_everystep=false, progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = ΔE
        rts[alg] = rt
        bytes[alg] = b
        allocs[alg] = memalloc
        nt[alg] = sol.destats.naccept
        nf[alg] = sol.destats.nf + sol.destats.nf2
    end
end
benchmark (generic function with 2 methods)

and change the timesteps. For adaptive integrators, we will use tolerances instead.

τs = 10 .^range(-4, -3, length=5)
ats = 10 .^range(-20, -14, length=5)
rts = 10 .^range(-20, -14, length=5)
5-element Array{Float64,1}:
 1.0e-20
 3.162277660168379e-19
 1.0e-17
 3.1622776601683793e-16
 1.0e-14

The tolerances for adaptive algoritms were chosen such that the energy error is close to the one given by the VelocityVerlet with the chosen τ.

Now, let us run the benchmarks

@progress "Variable dt" for (τ, at, rt) in zip(τs, ats, rts)
    runtime = Dict()
    ΔE = Dict()
    nt = Dict()
    nf = Dict()
    b = Dict()
    allocs = Dict()

    GC.gc()
    benchmark(ΔE, runtime, b, allocs, nt, nf, t, config(τ, at, rt))

    for (k,v) in ΔE
        push!(results, [string(k), runtime[k], at, rt, abs(v), nt[k], nf[k]])
    end
end

Let us now examine the results

results

35 rows × 7 columns

integratorruntimeτabstolEnergyErrortimestepsf_evals
StringFloat64Float64Float64Float64Int64Int64
1DPRKN12140.4871.0e-201.0e-200.001598579729177844
2Vern735.0461.0e-201.0e-200.0688734399841672
3DPRKN632.87941.0e-201.0e-200.0363926603545238
4Tsit518.28551.0e-201.0e-200.1444328121717
5VelocityVerlet583.3271.0e-201.0e-200.00332889350000700002
6McAte2588.8581.0e-201.0e-200.0001540773500001050002
7Vern9161.8571.0e-201.0e-200.099509211191191682
8DPRKN12141.3223.16228e-193.16228e-190.004434869853180346
9Vern735.03343.16228e-193.16228e-190.111018398441642
10DPRKN631.66583.16228e-193.16228e-190.0225537594044195
11Tsit518.45323.16228e-193.16228e-190.141024326421801
12VelocityVerlet326.4163.16228e-193.16228e-190.00405704196820393642
13McAte2331.7953.16228e-193.16228e-190.00058266196820590462
14Vern9162.4853.16228e-193.16228e-190.10493511173191618
15DPRKN12140.1641.0e-171.0e-170.01489219737177574
16Vern734.50651.0e-171.0e-170.0841265392540712
17DPRKN633.1741.0e-171.0e-170.0455213606946043
18Tsit518.48121.0e-171.0e-170.304543327421813
19VelocityVerlet183.2141.0e-171.0e-170.00762108110680221362
20McAte2184.5811.0e-171.0e-170.00368266110680332042
21Vern9162.9711.0e-171.0e-170.088116111225192226
22DPRKN12138.0893.16228e-163.16228e-160.007954249719177430
23Vern734.9453.16228e-163.16228e-160.112498395941112
24DPRKN631.86773.16228e-163.16228e-160.0207562589743922
25Tsit518.06513.16228e-163.16228e-160.108319324921465
26VelocityVerlet104.6513.16228e-163.16228e-160.0037017562240124482
27McAte2103.1373.16228e-163.16228e-160.0038468662240186722
28Vern9160.9113.16228e-163.16228e-160.097113310997188226
29DPRKN12142.7371.0e-141.0e-140.01563399849179842
30Vern734.73861.0e-141.0e-140.034158388440202
31DPRKN631.46331.0e-141.0e-140.100083587243768
32Tsit518.44251.0e-141.0e-140.314547325921783
33VelocityVerlet58.40341.0e-141.0e-140.003665623500070002
34McAte257.5481.0e-141.0e-140.0098037435000105002
35Vern9162.8731.0e-141.0e-140.12917311065190002

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

If we consider the number of timesteps instead, we obtain

@df results plot(:EnergyError, :timesteps, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of timesteps")

and we can also see how the muber of acceleration function calls affects the runtime

@df results plot(:f_evals, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Number of f calls", ylabel="Runtime (s)")